%% 体数据插值 volume data interpolation

%% 读入数据
clear; clc;
timerVal = tic; % 对整个脚本计时
[X, Y, Z, E1, nu1] = grdecl2fem('Grid_A.out','A_porps.txt');
[X2, Y2, Z2] = grdecl2fem('Grid_B.out');
elapsed_time_readdata = toc(timerVal);
fprintf('读入数据用时：%f s\n', elapsed_time_readdata)

%% 计算网格向量 grid vectors
xg = (X(1:end-1) + X(2:end))./2;
yg = (Y(1:end-1) + Y(2:end))./2;
zg = (Z(1:end-1) + Z(2:end))./2;

xg2 = (X2(1:end-1) + X2(2:end))./2;
yg2 = (Y2(1:end-1) + Y2(2:end))./2;
zg2 = (Z2(1:end-1) + Z2(2:end))./2;

% 单元属性值
Nxyz = [100 100 60];
E1 = reshape(E1, Nxyz);
nu1 = reshape(nu1, Nxyz);

%% 创建 griddedInterpolant 对象
tic % 对插值计时
F = griddedInterpolant({xg, yg, zg}, cat(4, E1, nu1), 'linear');
% 调用该对象完成插值
V2 = F({xg2, yg2, zg2});
elapsed_time = toc;
fprintf('插值用时：%f s\n', elapsed_time)


%% 将插值得到的杨氏模量、泊松比写入文件
% 格式化杨氏模量、泊松比 array
E2 = V2(:,:,:,1); nu2 = V2(:,:,:,2);
num_col = numel(E2)/6;
E2 = reshape(E2, 6, num_col);
nu2 = reshape(nu2, 6, num_col);

% 写入文件
fileID = fopen('B_props.txt', 'w');
fprintf(fileID, '\n\nYOUNGMOD\n');
fmt = [repmat(' 1*%.6f', 1, 6) '\n'];
fprintf(fileID, fmt, E2);
fprintf(fileID, '/\nPOISSONR\n');
fmt = [repmat(' 1*%.15f', 1, 6) '\n'];
fprintf(fileID, fmt, nu2);
fprintf(fileID, '/\n');
fclose(fileID);

%% 总用时
elapsed_time2 = toc(timerVal);
fprintf('总用时：%f s\n', elapsed_time2)